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TECHNICAL NOTE D-473 


AN ANALYTICAL INVESTIGATION OF THREE GENERAL METHODS 
OF CALCULATING CHEMICAL -EQUILIBRIUM COMPOSITIONS 

By Frank J. Zelesnik and Sanford Gordon 


SUMMARY 

The Brinkley, Huff, and White methods for chemical-equilibrium 
calculations were modified and extended in order to permit an analytical 
comparison. The extended forms of these methods permit condensed species 
as reaction products, include temperature as a variable in the iteration, 
and permit arbitrary estimates for the variables. 

It is analytically shown that the three extended methods can be 
placed in a form that is independent of components. In this form the 
Brinkley iteration is identical computationally to the White method, 
while the modified Huff method differs only slightly from these two. 

The convergence rates of the modified Brinkley and White methods are 
identicalj and, further, all three methods are guaranteed to converge 
and will ultimately converge quadrat ically. 

It is concluded that no one of the three methods offers any sig- 
nificant computational advantages over the other two. 


INTRODUCTION 

The determination of equilibrium compositions for systems of many 
constituents is generally difficult because the equations to be solved 
are not simultaneously linear. Since a direct solution is usually not 
feasible, some iterative technique must be used to obtain the solution. 

In the past few years there have appeared in the literature many articles 
dealing with chemical-equilibrium calculations of complex systems and 
describing various systematic iterative techniques (see list of refer- 
ences in ref. l). These articles present methods of solution applicable 
to specific chemical systems as well as general methods applicable to 
most chemical systems. Of the general methods available, those of 
Brinkley (ref. 2), Huff et al. (ref. 3), and White et al. (ref. 4) are 
perhaps the most widely used. Because of their wide use, these three 
methods were Investigated to determine whether any one of the three 
offers significant computational advantages. 
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The number of iterations and the amount of computation per iter- 
ation may be taken as the criteria of computational advantage. Since 
the amount of computation per iteration is essentially the same for the 
three methods, the only remaining criterion is the number of iterations. 
However, for any of the three methods investigated, the number of iter- 
ations is strongly dependent upon the initial estimates. Very poor 
estimates will generally require considerably more iterations than a 
good set of estimates. 

In order to make an analytical comparison, it is essential that 
all three methods begin with the same initial estimates. The Brinkley 
and White methods cannot start with the same estimates. The reason for 
this is that the Brinkley method requires the estimates to satisfy 
equilibrium conditions, whereas the White method requires them to satisfy 
mass balance. If a unique, real, and positive solution exists, these 
requirements are mutually exclusive except at the solution point. To 
permit identical estimates, these methods were modified to remove un- 
necessary restrictions on the estimates. 

In addition to these necessary modific itions, some other modifica- 
tions were made. These include simplifying the iteration equations of 
the Huff and Brinkley methods and treating condensed products in a man- 
ner different from that originally proposed for the Brinkley method 
(ref. 2) and White method (ref. 5). The Brinkley and White methods were 
also extended to permit using temperature as a variable. 

It will be shown that the three modified and extended methods are 
essentially equivalent computationally. 
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SYMBOLS 

total mass of reactant 
converged value of A 

formula numbers giving gram-atoms of element in specie 

gram-atoms of l^ 1 element per unit mass of mixture (eq. (3)) 

assigned value for gram-atoms of i^ 1 element per unit mass of 
reactant 

heat capacity per mole at constant pressure divided by R 
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formula numbers of components giving gram-atoms of i^ h element 
in j^ 1 component 

total free energy of mixture divided "by RT (eq. (7)) 

standard- state free energy per mole of specie = 

( H ?)j - T ( S T)j 

free energy per mole divided by RT for specie (eq. (8)) 
defined by eq. (68) 
defined by eq. (25) 
defined by eq. (48) 
enthalpy per mole of specie 
'fx/RT 

V RT 

defined by eq. (68) 

enthalpy per unit mass of reactant (eq. (15)) 

assigned enthalpy per unit mass of reactant 

number of different chemical elements 

number of gaseous reaction products 

total number of reaction products 

static pressure, atm 

assigned static pressure, atm 

partial pressure of specie, atm 

defined by eq. (44) 

universal gas constant 

number of reduced iteration equations 
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X j 

X J 

yJ 





A^ 



entropy per mole of specie in standard state 

entropy per unit mass of mixture divided by R (eq. (16)) 

entropy per mole of the specie divided by R (eq. (17)) 

assigned entropy per unit mass of mixture divided by R 

defined by eq. (27) 

absolute temperature 

assigned absolute temperature 

k th component of solution vector of modified Huff iteration 
equation 

k th component of solution vector of modified Brinkley iter- 
ation equations 

i^* 1 component 

total moles of gaseous products (eq. (4l)) 
moles of specie in a mixture 

moles of specie in a unit mass of mixture 

specie 

moles of specie in a mixture at equilibrium 

ith element 

activity of specie 
Krone cker delta 
defined by eq. (54) 
defined by eq. (68) 
defined by eq. (81a) 

moles of component in k^ 1 specie [see eq. (36)] 
i^b Lagrangian multiplier 
defined by eq. (34) 


fed 
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GENERAL DESCRIPTION OF ORIGINAL METHODS 

The three general methods all use an iterative technique to obtain 
equilibrium compositions. Initial estimates are made for the variables, 
and corrections to these estimates are obtained. The process is con- 
tinued until some arbitrarily selected criterion for convergence has 
been reached. 

Brinkley (ref. 2) was the first to treat the problem of the numer- 
ical solution of the nonlinear chemical-equilibrium equations for a gen- 
eral chemical system. In his calculation method the thermodynamic state 
of the system is specified by assigning the temperature and static pres- 
sure of the reaction products. The chemical-equilibrium relations are 
written in terms of components. Components are defined to be those in- 
dependent constituents of the mixture that can be used to express the 
overall composition. The number of components usually equals the number 
of different chemical elements in the mixture but under some rare cir- 
cumstances can be less (e.g. , if some elements appear in constant ratio 
in all the reaction products). The choice of components is not unique. 
The nonlinear set of equations defining the problem Is approximated by 
a linear set of correction equations obtained by a Taylor series ex- 
pansion of the nonlinear equations neglecting terms higher than first 
order. In Brinkley* s method any condensed reaction products that appear 
are treated as components. 

The method of Huff, Gordon, and Morrell was described in reference 3 
and subsequently presented in slightly modified form in reference 1. 

This method differs from Brinkley 1 s in the following respects: 

(1) The gaseous atoms are arbitrarily selected as components. 

(2) Condensed reaction products are not considered to be components. 

(3) The thermodynamic state of the system may be specified by 
assigning the pressure and either the temperature, the enthalpy, or the 
entropy. 

(4) The linear set of equations is obtained from a Taylor series 
expansion involving both logarithmic and linear variables. 

( 5 ) Corrections are applied to all constituents. 

The most recent of the three methods is that due to White, Johnson, 
and Dantzig (ref. 4). It is based on the criterion for equilibrium 
(cLF )rp p = 0. The thermodynamic sta^te of the system is specified by 

assigning the temperature and pressure. As originally presented, the 
method could only treat problems involving gaseous reaction products. 
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It was later extended to include condensed reaction products (ref. 5) 
by redefinition of one of the variables in the iterations equations. 


EQUATIONS FOR DETERMINING EQUILIBRIUM COMPOSITIONS 

The formation of any one of n chemical species from Z elements 
may be written in the form 

i a Z 1 = yj (l - 0 ~ n) (1) 

i=l J 

. , t 

where Z 1 is the symbol for the i^ chemical element, and Y^ is the 
symbol for the jth chemical species. Thus, a^ ^ represents the number 
of gram-atoms of element Z 1 in one formula weight of the chemical 
compound designated Y J; that is, the first subscript on indicates 

the chemical element and the second subscript gives the chemical com- 
pound. The range of any index such as k wil.i usually be associated 
with specific chemical species: 

1 < k £ Z gaseous elements 

Z + 1 < k < m gaseous compounds 

m + 1 S k < n condensed species (compounds and elements) 

where, for any particular problem, Z, m, and :i are to be considered 
as fixed integers giving the number of elements, the number of gaseous 
products, and the total number of all products, respectively. For the 
chemical elements (l £ j £ Z), 

a ij ~ ®ij 

where is the Kronecker delta. 

The overall composition of a mixture of :i chemical species may 
be expressed in terms of bj_, the number of gram-atoms of the i^ ele- 
ment per unit mass of mixture: 


n 

E 

j=i 


a i j x J 


(1 < i < l) 


(3) 


where xt gives the number of moles of the product in a unit mass 
J 

of the mixture. In some cases It is advantageous to alter the magnitudes 
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of the xl; and therefore, if equation (3) is multiplied by the total 
mass A, one obtains 

n 

Ab i = Yj a l.j x j (1-1 - l) (4) 

0=1 

where Xj is the total number of moles of the product in the mix- 
ture, and is given by 


w 


In the calculation of the equilibrium composition of a mixture , the 
number of gram-atoms of the i^h chemical element per unit mass of mix- 
ture is specified to be some assigned value b£. The condition for con- 
servation of mass then takes the form 


b? - b ± = Abi = 0 (1 < i < Z) (5) 

If y j are the values of xj that satisfy equation (5), then 

n 

Ab i = T, a ijyj = 0 (1 < i < Z) (6) 

J“1 


The condition for chemical equilibrium in a system at a tempera- 
ture T and a pressure P may be stated in two equivalent forms. The 
first says that the total free energy of the system divided by RT 


F(x J? T) = £ fjX. (7) 

0=1 


is a minimum at a constant temperature and pres sure , where 

+ In (l — j — n) 


f ..2n 

J RT 


0 


( 8 ) 


The values of Xj that minimize equation (7) at T,P are then the 
equilibrium values. In equation (7) and in all the following equations 
it is assumed that the standard state for gases is taken to be the 
ideal gas at 1 atmosphere, while for solids and liquids it is the pure 
solid or liquid at 1 atmosphere. It is assumed that the gas mixture 
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behaves as an ideal solution and that, if conlensed species occur, they 
will occur as the pure solids or liquids whos ; activities are independent 
of pressure; therefore, 


p. (1 < j £ m) 

1 (m < j £ n) 


(9) 


An alternative but equivalent way of stating the condition of chem- 
ical equilibrium is that the free-energy change across a reaction is 
zero: 


w 


Af j = 0 (1 < j < n) 

The free-energy change across the reactions (l) is 

l 

** j - f j - £ a u f i (i ~ j - n) 

i=l 

Nonequilibrium values of the composition variables will not minbnize F 
and will not make all Af j equal zero. 

The thermodynamic state of the system is specified by assigning the 
temperature and pressure. The static pressure of a mixture of gases is^ 
by Dalton T s law, 

m 

f-E'i (12) 

i=l 

The condition that at equilibrium the static pressure of the system is 
P Q is 

P = P Q (13) 

The specification of the thermodynamic state may be completed by re- 
quiring that, at equilibrium, 

T = T 0 (14) 

Equilibrium compositions are usually obtained for a specified pres- 
sure and temperature. For a number of problems, it is desired to obtain 
temperature and other system properties corresponding to an assigned 
pressure and some other assigned thermodynamic property such as enthalpy 
or entropy (e.g. , determining flame temperatures or conditions following 
an isentropic expansion). This is generally accomplished by calculating 


( 10 ) 

( 11 ) 
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data at several assigned temperatures followed by interpolation. How- 
ever, temperature can be determined directly during the iterative process 
for assigned values of enthalpy or entropy. This is accomplished by 
permitting temperature to be a variable and including an additional 
equation involving enthalpy or entropy. The choice of which technique 
is used - that is, fixed temperature for several points and then inter- 
polation, or temperature as a variable during iteration - is largely 
a matter of personal preference. 

The enthalpy of a mixture is 



where it is assumed that the enthalpy of the constituents is independent 
of the pressure and that reference values consistent with heats of 

) all n species. In like fashion the 

n 

= E s j x j (16) 

j=i 


In aj (l < j < n) (17) 

The temperature of the equilibrium mixture may now be specified either 
by equation (14) or by (18) or (19): 

£ = (. 18 ) 

S = s Q (19) 

where is an assigned enthalpy per unit mass, and s Q is an as- 

signed entropy per unit mass divided by the gas constant. 

The problem of chemical equilibrium is completely determined by 
either the set of conditions (5), (10), (13), and one of the three (14), 
(18), (19), or by the set of conditions (5), (13), one of the three 
(14), (18), (19), and the minimization of (7). Except for the simplest 
of cases, the set of equations selected cannot be solved in closed form 
and must be solved by iteration. Further, the equations are not all in 
the same composition variable; some are written in terms of the mole 
numbers Xj, others in terms of partial pressures p^, and still others 

in terms of both x^ and pj. The selection of the appropriate set of 


formation have been assigned t 
entropy of a mixture is 

As 


where 

Ui 

R 


s * - 
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equations, the choice of the iteration methoc to he used in their solu- 
tion, and the relation between Xj and pj that is used constitute the 
important differences in the three methods tc be discussed. 


MODIFIED METHOD OF HUFF ET A L- 

The iteration scheme of the Huff method is based on replacing 
the nonlinear set of equations by a set of l:.near correction equa- 
tions obtained by a Taylor series expansion neglecting terms higher 
than first order. The equations selected are (5) for mass balance, 

(10) for equilibrium, (13) for assigning the pressure, and either (18) 
or (19) for assigning the enthalpy or entropy. Although l + n + 2 
equations appear in this set, only n + 2 of these are independent, 
because with the use of (2) it is easily shown that the first 1 
equations of (10) are merely identities. 

If A is considered to be a constant, -.here are n + 2 equations 
in the m + n + 1 unknowns pj, Xj, and T. However, if A is treated 

as a variable, there are then m + n + 2 unlmowns, and it is possible 
to adjoin to the previous set of equations n additional equations re- 
lating Pi to x v In particular, the set selected is 
J J 

Pj = Xj (l - j 5 n) (20) 


These m equations may now be used to eliminate Pj from the previous 
set of equations, resulting in n + 2 independent nonlinear algebraic 
equations in the n + 2 variables Xj, A, and T. If these are ex- 
panded in a Taylor series about an estimate for the solution using 
In Xj(j £ m), xj(m < j < n), In A, and In T as the variables, one ob- 
tains the following equations linear in the correction variables 
A In Xj, Axj, A In A, A In T from the set o? nonlinear equations 

(5), (10), (13), and (18): 


= A Abi 



m 

n 

n 

= E a ij x j A ^ x j 

+ E a u Ax j • 

- E 

j=i 

j^m+1 

j=i 


"Af j = A In Xj - £ a kj A 111 x k " ( h j " E a kj h kV 

k=l \ k=l / 


( 21 ) 

In T (1 < j < m) 


(22a) 
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If y v V and. T Q are the equilibrium values of x-, A, and T, then 
J J 

A In xj = In yj - In xj (l £ j £ m) 

Axj = yj - Xj (m < j < n) 

A In A = In A 0 - In A 
A In T = In T q - In T 

For xj, A, and T sufficiently close to the equilibrium values, these 
linear equations will give the equilibrium values. In general, however, 
the values obtained will only be better estimates for the equilibrium 
values than xj, A, and T. 
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If the same procedure is followed for (19) and if equation (23) is 
then added to the result , the following equation is obtained: 

A(s q - s) + AP = A As + AP 


m 


= £ SjXj Aklj, £ si AXJ - £ sJ Xj A to A 
j=l j=m+l j=l 


where 


+ ^ CjXj A In T (26) 
j=l 


^(s°). 

- ft J " 111 x j (l < j < m) 


* 

S j " < 


R 


(27) 


(m < j < n) 


Equations (2l), (22) ; (23)^ and either (24) or (26) can readily be 
reduced to the following Z + (n - m) + 2 linear equations by using 
(22a) to eliminate A In and (25) to eliminate Af?: 


r ik a k + X* a ik ^k + ^ a ik x k(~^ 
k=l k=m+l k=l 


k k=l 


/ m Z \ m 

+ f X^ a ik h k x k ’ X* r ik h k ) & In T = A Ab.^ + ^ a ik x k f ' 

‘ “ k=l 


k=l 


* 

- Y r ik f k (1 a i < 1) 


(28) 


k=l 


l 

£ 

k=l 


a kj °k + ' h < 


? - £ a kJ h°\ a k . f J 

k=l / k=l 


(® < j - n) 




( 29 ) 
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A In T 
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Z / m 

EE 

k=l V .1=1 


E 

j=l 



^ • E f 


j=l 

n 

n 

E 

hg Ax k + 2 h ! 

k=m+l 

j=l 

m 

l / m 


Z / m 


a kj x j\ f k 


k=l\j=l 


C j X j 


a kj h j x j| h k 


A In T 


m 


l / m 


A ^ + E h j f .j x j E(E a k.j h j x j) f k 


,i=l 


k=l\,i=l 


I / m 

z(z 

k=l\.j=l 


a kj ! 


| °k + E s J 

+ E 3 J x 5 ( " A 


/ 

,j=m+l 

j=l 


n 

m 

l / m 

+ 

E 

c j x j + E 

- E E 



rH 

1! 

•r-y 

k=l\,i=i 


m 




= A As + AP + £ Bjf jxj 

j=l 


Z / m 


-EE 

k=l\j=l 


A In T 


* \ 

a kj s j x j] I k 


//here 


m 

j==l 


Hj a kj x j 


(l £ i <1, 1 £ k < l) 


J k 


= A In x-^ ( 1 £ k £ m) 


(30) 


(31) 


(32) 


(33) 

(34) 


The set of Z + (n - m) + 2 linear equations (28) to (30) and 
either (31) or (32), which can he solved to give directly the corrections 
to the estimates for yj (l < j < l, m < j < n), A 0 , and T 0 , is pre- 

sented in figure 1. These corrections can also he used to obtain the 
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corrections for the other gaseous species from equation (22a) rewritten 
using (34) and (25): 

= + a kj(fk + a k) + “ V a kjh§JA In T (l — j — m ) 

k=l \ k=l / 

(35) 

The equations for (l — j — l) in (35) are identities for the gaseous 
atoms hut are included for the purpose of later comparisons. 

When the temperature is assigned directly hy (14), temperature is 
no longer a variable, and the appropriate set of iteration equations is 
(28) to (30) with the A In T terms deleted. The A In T term is also 
deleted from (35). 


MODIFIED METHOD OF BRIMLEY 

In the previous section the equations for conservation of mass and 
for chemical equilibrium were written in terms of the gaseous atoms. The 
corresponding relations in the Brinkley method are written in terms of 
components. A set of iteration equations in terms of components is pre- 
sented first. These equations are then converted to an equivalent set 
of equations in terms of gaseous atoms. 


Iteration Equations in Terms of Components 

If the component is designated as then all of the chemical 

species in an l -element system can usually be written as a linear com- 
bination of l components as follows: 

V Vj k wj = Y k (1 < k < n) (36) 

j=l 

where gives the number of moles of the component in the k^* 1 

species^ may be expressed in terms of the elements as 

= £ c^z 1 (1 < j < 1) (37) 

i=l 

where c j gives the number of gram-atoms of the i^* 1 element in the 
j'k* 1 component. Substituting (37) into (36) end comparing coefficients 
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of with those in (l), 

T. c ij v jk = a ik (l ^ i — 1 £ k £ n) (38) 

j=l 

To simplify the discussion, the last l chemical species 
(n-7,+1— k — n) will be selected as components; that is, the j^ 1 
component is the (j + n - l)^ species, or 

WJ = yj+n-Z (l < j < l) (39) 

This selection includes the condensed species as components as required 
by the Brinkley method. This particular choice of components will in 
no way restrict the discussion, since it will subsequently be shown that 
the choice of components does not affect the corrections obtained during 
the iteration procedure. Because of the particular choice of components 
(39), 


c ij — a i, j+n-Z U — i — n 1 — J — ^ ) (40) 

The vjk can be obtained from (38) if the matrix Cj_j is nonsingular 
and possesses an inverse c^. The requirement that c-^j be non- 
singular restricts the possible choices of components. In the previous 
section, because the atoms were used as components, c ij “ s ij and thus 
was always nonsingular. 

In the previous section, it was found convenient to treat A as 
variable in order to introduce a simple relation between the moles and 
partial pressures (20). In this section, A is regarded as an arbitrary 
constant A r , and a new variable x is introduced. The variable x is 
equal to the total moles of gaseous products: 

m 

x = Xj (41) 

j=l 

The partial pressures are then related to the moles by 


P 


j 



(1 £ j < m) 


(42) 


The equations for conservation of species are readily obtained. 
Substituting (38) into (4), multiplying the resulting equation by c“^, 
and summing on i from 1 to l give 


A,< ik = ^ 




(1 < k < Z) 


(43) 
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where 


3-k 


t 

i =1 


c ki^i 


The conservation equations corresponding to (5) are 

q£ - q. k = Aq k = o (i-k-i) 


(44) 


(45) 


The condition for equilibrium is again taken to he (l0)j however , 
the free-energy change across the reactions (36) with the choice of com- 
ponents (39) is 


l 

^k = *k - E V jk f j+n -l (1 < k < n) (46) 

j=l 


Since the last l species were selected as components, 

v jk = Sj+n-Z,k U-tl,n-J<k<n) (47) 

and thus the last l equations of (46) are identically zero. 

In his original presentation, Brinkley eliminated the variable x 
from his equations by imposing the condition 


l 

A' E 


k=l 



1 




For comparison with the other methods, the variable x is retained. 
The use of equation (42) ensures that (13) is satisfied for any choice 
of xj. Equation (13) is replaced by the relation 


G = 



- 1 


(48) 


If equation (41) is satisfied (i. e. , if during the iteration the esti- 
mate for x is obtained from the sum of xj(l £ j £m)), then 

G = 0 (49) 

If, however, an independent estimate is made for x, then G will be 
zero only when the iteration converges. 
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The equations that must be satisfied by the system in equilibrium 
are (45), (10), (49), and either (14), (18), or (19). Using (42) in (10), 
expanding the resulting equations together with (45), (49), and (18) in a 
Taylor series a;s in the previous section, and treating (xj/x)(l < j < m), 
Xj(m < j 5 n), x, and T as variables give 

A' (q.i - q ± ) = A' Aq^ 


m 


m n 

= £ v ik x k^k + £ v ik ^k + 2 v ik x k ^ (l < i < l) 

k=l k=m+l k=l 


-"k “ 5k - £ v jk £ 


7+m-n 

E 

J=l 


jk^j+n-Z - h k 


- i f 

j=l ) 


(50) 

(l < k < m) 
(51a) 


1 + m-n 

L k Zy v jk^j+n-Z 

J-l 


-Afv = - 


" ^ h k " ]C v jk h j+n-Z^ 


o \ Zff / i <r \ 

1 — (m < k E n) 


m m 

X " Yj X J = Yj Xk ^ k 

J-l k=l 


(51b) 

(52) 


A'(h Q - h) a A r Ah 


m n m — n 

7! h k x k^k + h k ^k + h k x k C k x k V 

k=l k=m+l k=l x k=l 


where 


?k * 


(y k /y) - ( x k/ x ) 


x k/x 

^k = yk - x k 

Ax = y - x 


AT = T 0 - T 


(l E k E m) 
(m < k < n) 


y ( 54 ) 
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Treating (19) in similar fashion and adding (l + In ») times (52) to 
the result give 




A T As + II + In 


( 1 + to T) /_ 


m \ 

- E x j) 

m 

= Y s ^ Xk ^ + 

n 

E 

C_i. 

II 

H 

k=l 

k=m+l 

/ m 

p \ — 

n 

( E s A 

\k=l 

- ; m *£ 

X J X 

- E ( 

k=l 




AT 


Assuming that equation (4l) is used to obta:.n an estimate for x, the 
left side of (52) and the corresponding tern on the left side of (55) 
are zero. These terms will therefore he dropped from further considera- 
tion. However, it is possible to treat x in a manner identical to the 
other variables and to alter the estimate for x as called for by the 
iteration equations, so that (4l) need not be satisfied until equilibrium 
is reached. Equations (50), (5l), (52), and either (53) or (55) can be 
reduced to Z + 2 independent equations by using (51a) to eliminate 
from the other equations and (46) to eliminate Af^ as follows: 


Z+m-n / m \ n 

X [£ v ik v jk x kUj+n-z + X V;k 

j=l \k=l J k=m+l 


Axv + 


m 

Y v ik x k 


k=l 


Ax 

x 


m 


r.i 


S v ik h k x k “ ^(^2 v ik v jk x kj h j+n-Z 
k=l j=l\k-l / 


AT 

T 


Z-Hn-n 

E 

j-1 


- A'Aq i + X v ik f k x k -tit v ik v j k x k^ ^ j +n - Z 

j : : 1 \k” 1 / 


k=l 


Y v jk£j+n-Z + ( h k 


Z 

E 

j-i 


' X v jk h j+n-l)^F “ f k - X V jk i j+n -l 

/ j-1 

(m < k E n) (57) 


(1 < i < l) (56) 

l 


7+ra-n./ m \ rm ‘ / m \ 

( X V jk x k)^+n-Z + X h k X k ~ X(X v jk x k] h j+n-Z 
j-1 \k=l / Lk=l j“l\k=l / 


AT 

T 



( 58 ) 
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D 

H 

i 

q 


fm-n/ m \ n m — 

2(2 h k v jk x k)£j+n-Z + ^ h k ^k + 2 h k x k £" 

3=1 \k=l ) k=m+l k=l 

[ • n m Z / m \ -i 

52 °k X k + 2 h k h k x k " £(52 h k v jk x kVj+n-Z 
k=l k=l j=l\k=l / J 


AT 

T 


m 

= A* Ah. + S S 


E 

k=l 


Z / m 


- E E * 


jk x k| f j+n-7 


(59) 


j=l\k=l 


«n-n/ m v n m / p \ — 

E (E ^jaVj«-z ♦ E + E v* - ;i» f )f 

j=l \k=l J k=m+l k=l ' 

"n m Z/m \ 

+ E °k X k + 52 S k h k x k “ 2(52 s k v jk x k] h j+n-Z 
k=l k=l j=l\k=l / 


m 


Z/m 


= A'As + £ s k f k x k - 52 ( 22 S k v jk : 


AT 

T 


" vX kKj+n-Z 


k=l 


j=l\k=l 


(60) 

Equations (57) are identically satisfied and are included only for 
later use. The 7+2 linear equations (56), (58), and either (59) or 
(60) can he solved to give directly the corrections to y^n - l <k^n) 

and T. The corrections to the estimates for y^(l £ k £ n - l) are 
obtained from (51a). Equation (51a) can be written in a form similar 
to (35) as follows: 


£k 


+ In 



x 


7-Rn-n l 

- f k + 52 V jk£j+n-Z + 52 V jk f j+n-z 

j=l j=l 



(1 < k < m) 


(61) 


The equations for the gaseous components (n - l < k < m) in (6l) are 
identities. 


When the temperature is assigned directly by (14), temperature is no 
longer a variable, and the appropriate set of iteration equations is (56) 
to (58) with the AT/ T terms deleted. The AT/T term is also deleted 
from (6l). 
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Iteration Equations in Terms of Atoms 

The set of l + 2 independent equations (56), (58), (59), plus the 
n - m identities (57), or the same set of equations with (59) replaced 
by (60), can be converted to l -f (n - m) + 2 independent linear equa- 
tions very similar in form to equations (28) to (32) used in the Huff 
method as follows: 

First, the may be obtained from (38) as 


v jk ~ ^ c ji a ik 
i=l 


(62) 


Next, equations (62) and (44) are used in (58) to (60) and the equations 
resulting from (56) are also multiplied by c and summed on i. The 

resulting system of equations is: 


7 

r 

k=l 


n 


Y r ik£k + Y &ik ^ + E a ik A k 


m / m 

a-iOiv — + V a-rvhvx 


k=m+l 


k=l 


- + E 


,, 0 .* 

ik u k x k 


E 

k=l 


h o\AT 
r ik h k|7jT 


/ 

A’Aki + Y, a ik x k( f 


f k + In 


k=l 


V 




l 

k=l 


r ik^k 


Y &k 0 ^ k + " E ak J hk ] T 8 ~ E a k, 

k=l \ k=l / k=l 


(1 < i < 1) 

(63) 


(m < j — n) (64) 


l / m \ rm I 'm \ 1 A m 

]£[]£ ak J x jV k + E h S x j ■ E(E t 

k=l\j=l j Lj»l k=l\j=l / . 


0=1 


P Z / m 

■ e(e 

k=l\j=l 


= Y f !h +xlI1 f ■ E(E ^koVk 


(65) 
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Z / m 


EE + E 


m — 

h 2 ^ + E h h f 


k=l\j=l 


k=m+l 


j=l 


n m Z / m 

E c j x j + E ' E(E 

•j=l j~l k=l\j=l 


,o \r«|AT 
av i I hi. 


a kj j jl k| t 


m / p \ i/m v _ 

= A' Ah + ^ hjXj/fj + In -S\ E a kj h j x jV: 

j=l V 7 k=l\j=l / 


1 

l / m 


(66) 


E ( t + E ^ t) f 


k=l\j=l 


k=m+l 


^k=l 


o \ Ax 

X I X 


“ n 111 2/m \ a m 

E c j x j + E s K x j - E(E f 

J=l o=l k=l\j=l / _ 


ID 


E s l x j ( f J 

o=i 


P \ Z / m 

f) ■ £5 


a kj s j x j] f k 


(67) 


k=l\j=l 


where 




cTi f . 


E "Jk ^J+n 

0=1 


__ Z+m-n 

^k = E C 0 ^ ^j+n-Z 

0=1 


^k ^ c jk *\i+n-l 
J =1 

2 p 2+m-n 

, = V c:i f* + In V cTi 1 

-l L, Jk j x Z-< J k 


0=1 


j=l 


Hl<k<Z) (68) 


J 


In similar fashion, equation (6l) becomes 


P Z 

r o „* 


+ 1X1 x “ f j + E ^j^k + ^ + h j " E T 


k=l 


Z 

E 

k=l 


(l < j < m) 


( 69 ) 
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The set of equations (63), (64), (65), and (66) or (67) represents 
a more useful form of the iteration equations for automatic computation 
than the set (56), (58), and (59) or (60). Although the latter contain 
n - m fever equations if condensed products are present, it is necessary 
to obtain a new set of for each nev choice of components, and the 

calculation can be rather lengthy when many species are under considera- 
tion. However, using the former set it is only necessary to calculate 
an inverse matrix for each new choice of components and from it the f k 

and h£. In a subsequent section it will be snown that it is possible 
to obtain a set of iteration equations that do not contain f^. and h£, 
and hence it is unnecessary to obtain the matrix c?j. 

MODIFIED METHOD OF WHITE 3T A L. 

The methods of Huff and Brinkley are Newton-Haphson iterations 
where the conditions for chemical equilibrium are given in terms of free- 
energy changes across reactions (equilibrium constants). The White 
method uses the alternative but equivalent method of specifying equilib- 
rium conditions as the minimum of the total fr^e energy of the mixture. 

With (42), equation (7) becomes 

n _ p 

F(*i,T) = £ f*x ± + x Ln ~ (70) 

i-1 X 


If (yj_jT 0 ) represents a neighboring point to (xj_,T) with y^ satisfying 

equation (6) with A = A*, then the Taylor serias approximation for the 
free energy at (yj_^T 0 ), expanded about (xp,T) and neglecting terms higher 
than second order, is 

m / P \ 1 

«(7i,T 0 ) = F( Xi ,T) + £ (f * + In + Y f*Ax ± 

i=l i=n+l 


E h?Xi f + 1 E En* - iW ** 


- i=l 


i=l k=l 


h?Ax,f + i£( Z ,J-cg)x k (f) 2 (71) 


k=l 


k=l 
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because 

BF + T i 

i [ft (m < i < n) 

BF _ y-' ^i x i 

§T " ' A* ~ T~ 

1=1 

and 

( 1 5 i 1 m ; 1 £ k < m) 


(m < i n, m < k 5: n) 


(1 < i < n) 

y ( 2h ° - c °>*i 
i=l 

The function Q is to he minimized as a ftu 
to the restriction 

A*Ab ± - £ aijfyj - x j ) = 0 (l 

j=l 

obtained by the subtraction of equations (6) and 
it is necessary to minimize the function 

u - « + i i a u<yj 

i“l L J=1 

where the are Lagrangian multipliers. This 

lowing conditions be satisfied: 

tT~ = ■ V *i a ik = 0 (l<k<n) (74) 

^ L> 

i=l 


iction 


of y^ subject 


< i < i) 


(72) 


( 4 ) with A = A 1 . Thus , 
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(73) 


requires that the fol- 
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= J x i X 




d 2 F 

Sx^T 
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From equation (7l), 
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+ /A - Xl - hfi ^ 


^ f 1 £ + In ^ + 

X X \ X^ X 
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1 f* _ h° — 

rp 

Thus, equation (74) becomes 


(l S k £ m) 
(m < k £ n) 


K x 


X - h o ^ . 

x k x/ k T 


Jti a i^ = 0 (l £ k £ m) (75a) 


i=l 


■p* v.O 
f k " h k rp 


Y - i, *i a ik " 0 (in < k < n) 

i=l 


(75b) 


Equations similar in form to the correction equations in the Huff 
method [(28), (29), (30), and (3l) or (32)] aid the Brinkley method 
[(63), (64), (65), and (66) or (67)] may now be obtained as follows. 
Rewriting, equation (72) in the form 

m n 

^ a jk^k “ x k) + ^ a jk ^k “ A T Abj 

k=l k=m+l 

and eliminating y^ (l £ k £ m) with equation (75a) give 


E + E 


n m 

a ik ^k ^ vX " 


k=l 


k=m+l 
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l ik A k Y + 
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E 

k= L 


, o 

a ik^k x k t 


m I P \ 

= A'Abp + £ a ik x k f k + ^ Y (l ~ 1 - Z) 

k=l V 


(76) 


Solving equation (75a) for y^. and summing over k(l < k < m), 
l / m 

t 

k=] 


Z / m \ m Am i , p 

E(E vjW + E “&■* f - 1 f k\ + X in ^ 

k=l\j=i ) k=i k=:_ 


( 77 ) 
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The remaining iteration equation can be obtained from either (18) 
or (l9). Expanding (18) in a Taylor series about (xj,T) gives 


A’Ah = J hj Ax k + ^ C°x k 


,o__ AT 
T 


k=l k=l 

Eliminating y k (l < k < m) from (78) using (75a) gives 


l / m 

E(E 

k=l\j=l 


a kj h j x j K + 
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k=m+l 


h k ^k + 
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h£x 


Ax 
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(78) 
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+ Z h k h k x k + Z C k x k]^ = A,Ah + Z h k x k f k + m -£) (79) 

\k=l k=l / k=l 


I 


Jr 

Proceeding in an identical fashion with (19) and adding In S times 
(77) to the result give x 


l / m 


E(Ev;4v E 

k=l\j=l / k=m+l \k=l / 


m 


n \ m , -r, \ 

E + E ■ A,As * E 3 ^( f i + ln f) 

k=l / k=l ' ' 


(80) 


,k=l 


Equations (76), (75b), (77), and either (79) or (80) are the 
iteration equations of the extended White method in the variables 
w k (l — k £ l), Ax k (m < k E n), Ax/x, and AT/T. The new values for 
khe Xj (l E j E m) are obtained from (75a) with the use of the solution 

to the preceding set of equations. Equation (75a) can be written in a 
form similar to (35) and (6l) as follows: 


Aj + ln -=■ = -ft + 


x 


+ hi — 


J ■ Z a kj*k “j T 

k=l 


(l < j < m) (81) 


where 


Ax^ 


J Xj 


Ax 


X 


(81a) 



COMPARISON OP MODIFIED METHODS 


Differences in Presentation from Original Reports 

Before making a comparison of the three extended methods, it is 
appropriate to point out how the methods presented here differ from the 
original presentation. The Huff et al. method as presented here is sub- 
stantially the same as that in reference 1, differing primarily in 
notation, a modified entropy equation, and the use of alternative ex- 
pressions for the elements of the last two columns of the augmented 
matrix. 

Brinkley expressed his correction equations in terms of £k* In 
this paper the correction equations are in terms of £k* Both sets of 
equations lead to the identical corrections. The use of permits 

a direct comparison with the other methods as well as offering some 
computational advantages, as was pointed out in a previous section. 

Brinkley’s original method was an iteration for assigned tempera- 
ture and pressure where composition estimates were made only for the 
components, the remaining composition variables being determined from 
the equilibrium constants. As presented here, the method may be used 
for variable -temperature iteration, and the estimates for all the com- 
position variables may be made independently. If, by choice, independ- 
ent estimates are made only for components, then the free-energy terms 
on the right side of the iteration equations vanish identically, and 
only the mass-balance and enthalpy or entropy errors remain. However, 
this latter procedure may lead to unnecessary difficulties, as for ex- 
ample when some of the components are not major species. In this case 
small errors in the trace components are magnified to such proportions 
that convergence may be slow if not impossible. 

The extended White method permits variable -temperature iteration, 
inclusion of condensed reaction products, and the use of estimates that 
need not satisfy mass -balance restriction. When the temperature of the 
system is assigned directly (14), only gaseous products are considered, 
and the estimates satisfy the mass-balance equations, then these equa- 
tions reduce to those of reference 4. For the case of variable tem- 
perature, the minimization procedure has been combined with a Newton- 
Raphson iteration for temperature. 


Comparison of Modified Iteration Equations 

The equations for the three modified methods are presented in fig- 
ure 1. In this comparison the set of equations for the variables 
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will be used for the Brinkley iteration, since these equations charac- 
terize the Brinkley iteration as adequately as the equations for the 
variables. A comparison of the corresponding Huff and Brinkley equa- 
tions discloses a great similarity; the coefficient matrices are identi- 
cal, except for the last two columns. The second-last columns differ 
only in that in the Huff equations the summation extends from 1 to n, 
while in the Brinkley equations the summation is only to m. This is 
directly attributable to the different roles played by the variables A 
and x. The former is associated with all the species, while the latter 
refers specifically to the gaseous reaction products. The last columns 
of the coefficient matrices are both formed in part from linear combina- 
tions of the first l columns, differing only in the fact that constants 
of combination differ. The constants of combination are h£ in one 
case and h£ in the other. The right sides also differ by the use of 
different constants of combination in forming linear combinations of 
columns of the coefficient matrix. 


» 


FINAL MODIFICATIONS 

Derivation of Component-Independent Forms 
of Huff and Brinkley Equations 

The elimination of the linear combination terms from the Huff 
and Brinkley equations in figure 1 would make these equations in- 
dependent of components. To eliminate these terms, it is necessary 
to establish a relation between the solution vector of a set of 
equations with linear combination terms and the solution vector of 
a new set of linear equations with these linear combination terms 
removed. 


To accomplish this objective, consider two sets of r linear 
equations: 


Wi j u j _ di 




0=1 

r 


> (l<i<r) 


E Mi j v j " 61 

j=l 


(82a) 


(82b) 
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where the relations between the matrices and N-y and between 


the vectors dj_ and e^ are 


M ij - ~ 6 jr N ipPp 

P=1 


i 


(83a) 
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= a. 


Z< N ip Y p 

P=l 


(83b) 


The constants (Bp and Tp are defined to le zero for p greater than 

l and nonzero for (l £ p £ l). The equations (82b) are analogous to 
the Huff or Brinkley equations of figure 1, while (82a) represents these 
equations with the linear combination terms deleted from the last column 
of the coefficient matrix and the right side. 

To obtain the relation between u^- and v^ ; both (82a) and (82b) 
are multiplied by the inverse matrix N~;j- end summed on i. Using the 
identities (83a) and (83b ), one obtains 


r 

u k “ ^ ^ki^i 
i=l 

r 

y k - v r^k = Yj N ki d i " r k 
i=l 

and therefore 

Y k = u k + v rPk " r k (l < k < r) (84) 

From (84) and the definition of the constants of combination there 
follows 




vk = 


u k + Y r^k " "^k 
u k 


(i e k e i) 

(l <k Sr) 


(85) 


Equations (85) is the desired relation. It shows that the linear com- 
bination terms affect only the first l components of the solution 
vector, leaving unaltered the remaining components of the solution vec- 
tor. The following table indicates the correspondence between the 
variables in (85) and the variables of the modified Huff and Brinkley 
equations of figure 1: 


QT 
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Modified 

Huff 

Modified 

Brinkley 

r-s 

VI 

X 

VI 

rU 

£ 

a k 

Ik 



M 1 

V r 

A In T 

T 

h. (1 < k < l) 

u 0 

h k 

ro 

h k 

r k (i < k < i) 

1 k 

? k 


With the correspondence of the table and the use of u^. - with 
the Brinkley variables, (85) for (l 1 k ^ l) gives the two transformation 
equations 

° k + f ^ = u. k + A in T (1 5 k Si) (86) 

Ik + = u k + Y (l - k < l) (87) 

Eliminating o k from the modified Huff equations (28),_(29), (30), (3l), 
and (32) with the transformation (86) and eliminating from the 

modified Brinkley equations (63), (64), (65), (66), and (67) with the 
transformation (87) give a new set of iteration equations for these two 
iteration methods. These are presented in figure 2, together with the 
modified White iteration equations from figure 1. 


All three sets of iteration equations in figure 2 now give directly 
the corrections to the condensed species, temperature and either x or 
A; however, none of the iteration equations give the corrections to the 
gaseous species directly. The corrections to the gaseous species in 
the modified White iteration are obtained from (8l). Corrections to 
the gaseous species for the Huff and Brinkley equations of figure 1 are 
obtained from (35) and (69), respectively. These two equations can he 
transformed to correspond to the variables of the Huff and Brinkley 
equations of figure 2 by substituting (86) into (35) and (87) into (69): 


°J = - f 1 + E + h° A In T 

k=l 


Pfh * 


Ej + * f - - f 1 '* £ + h J T 
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> 

k=l 
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(1 < j <m) 

(88) 

s 

VI 

•<-3 

VI 

1 — 1 

(89) 
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Comparison of Final Equations 

An examination of figure 2 discloses that in this form all three 
sets of iteration equations are independent oF components. Further, in 
this form the Brinkley and White equations are computationally identical, 
while the Huff equations differ only slightly from these two. The Huff 
equations have additional terms in the second-last column of the coef- 
ficient matrix, while the other two have additional terms on the right 
side. When only gaseous products are considered, the three coefficient 
matrices are identical except for one term in the entropy row. 


OBTAINING AND APPLYING CORRECTIONS 

A distinction should be made between the solution of the equations 
to obtain corrections and the application of these corrections to ob- 
tain new estimates. 


Obtaining Corrections 


As was pointed out previously, all three modified methods give the 
corrections to the condensed species, temperature and either x or A 
directly. The corrections to the gaseous species are obtained from the 
auxiliary relations (8l), (88), or (89)- For the Huff iteration, (88) 
gives the correction variable a-, which may be regarded as either 
Ax • J 

A In Xi or — — and applied accordingly. However, (8l) and (89) give 

J x j 

the variables Aj and respectively. 


Ax, 


The corrections A In Xj or may bs related directly to Aj 

or if use is made of the approximate relation A In z = — : 

J 2. 


Ax, 

A In x, = — si = y 
x, A 




. Ax 
A 1 + — 
J x 


t. +$£ 

X 


(modified White) 
(modified Brinkley) 


(90) 


Since the Brinkley and White iteration equations in figure 2 and also 
the corresponding auxiliary relations (89) and (8l) are identical, it 
follows that the modified Brinkley and White equations provide identical 
corrections. These corrections will in general differ from those given 
by the Huff iteration. 
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Applying Corrections 

Once the corrections have been obtained from one of the three iter- 
ation procedures, there still remains the question of how to use these 
corrections. Some of the possible variations are 

(1) Using the corrections for the gaseous species and temperature 

Ax k 

linearly or logarithmically, for example, or A In x k 

x k 

(2) Using the entire correction or some fraction of the correction 

(3) Using some or all of the composition corrections: 

(a) Applying only l of the composition corrections and de- 
termining the remaining n - l composition variables from 
the equilibrium constants 

(b) Applying n - l of the composition corrections, deter- 
mining the remaining l composition variables from the 
mass -balance relations 

(c) Applying all n composition corrections 

(i) Using some combination of the above three. 

In the original papers Huff et al. advocate the use of logarithmic 
corrections while Brinkley and White use linear corrections. Further, 
Brinkley makes use of variation (3a), and White (3b), while Huff uses 
(3c). In addition, White suggests the use of fractional corrections 
and reference 1 indicates that fractional corrections may also be used 
with the Huff method. 


CONVERGENCE OF THE THREE MODIFIED ITERATION METHODS 

Because of the experience factors involved in the use of correc- 
tions, no completely analytical discussion can be given to the problem 
of convergence rate. It is possible, however, to draw some conclusions 
concerning the relative convergence rates of the three modified methods 
if it is assumed that all three methods use the corrections in an equiv- 
alent manner. The modified Brinkley and Huff methods axe Newton-Raphson 
iterations, which are known ultimately to converge quadratically (ref. 

6). Since the modified Brinkley and White methods are identical, all 
three methods must ultimately converge quadratically. Further, Gleyzal 
(ref. 7) has demonstrated that the Newton-Raphson iterations are guar- 
anteed to converge if an appropriate use of the corrections is made; 
thus, all three methods are mathematically guaranteed to converge. 
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The remaining small differences between the modified Huff and the 
other two modified methods is due to the use cf the variable A in one 
case and the variable x in the other. It wculd therefore seem reason- 
able to expect that the rate of convergence would be about the same for 
all three. This is somewhat substantiated by the fact that, for the 
particular problem used as an example in reference 4, both the Huff and 
White methods required six iterations to converge to the same degree of 
accuracy. This was so even though logarithmic corrections were used in 
the Huff iteration and fractional linear corrections were used in the 
White iteration. While one example cannot be considered conclusive, it 
is indicative of the fact that none of the methods can be expected to 
offer a marked advantage. 

Experience has shown that by the proper use of the corrections it 
is possible to accelerate the convergence. Thus, if one is in a situa- 
tion where the errors are primarily in equilibrium, the use of loga- 
rithmic corrections will reduce the error most rapidly. This is true 
because the equilibrium equations are linear in the logarithmic variables. 
However, where the error is primarily in the mass balance, the linear 
corrections are best because the mass -balance equations are linear in 
the variables. Linear corrections permit negative values of xj, whereas 
logarithmic corrections do not. While negative Xj for condensed re- 
action products are significant since they indicate that the particular 
condensed product is absent, for gaseous products the xj should be 
greater than or equal to zero. 

Several years of experience (ref. l) have shown that the use of 
logarithmic corrections for the gaseous products and linear corrections 
for condensed products permits rapid convergence for a wide variety of 
problems even with poor first estimates. 


SUMMARY OF RESULTS 

Three widely used methods of equilibrium calculations were analyt- 
ically compared. To permit the comparison, the methods were modified 
and extended. The Brinkley method was extendel to allow variable tem- 
peratures and to permit estimates that do not satisfy the equilibrium 
constants. The White method was extended to accommodate condensed 
species and variable temperatures and also to permit estimates that do 
not satisfy mass -balance requirements. The comparison of the methods 
showed the following results: 

1. The modified methods are computationally equivalent. 

2. The modified Brinkley and White iterations give identical 
corrections. 
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3. Alternative sets of iteration equations for the Brinkley and 
Huff iterations were obtained that offer computational advantages over 
the original set. 

4. Neither the concept of components nor the particular set of 
components selected plays a role in obtaining corrections to the 
estimates. 

5. All three methods are guaranteed to converge and will ultimately 
converge quadratic ally. 

6. No one of the methods offers any significant advantage over the 
other two. 

Lewis Research Center 

National Aeronautics and Space Administration 
Cleveland, Ohio, July 7, 1960 
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Figure 2. - Comparison, of modified iteration equations with linear combination terms deleted (W, White} B, Brinkley} 
H, Huff). 
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